//#include <string.h>

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"

#include <vector>
//#include <deque> 

#include "Data.h"

using namespace std;


class ICs { 
	public:
	
	typedef unsigned int uint;
	typedef long unsigned int luint;
	typedef double real;
	
	Data &data;
	
	class Grid{
		public:
			int xn, yn, zn;
			double xs, ys, zs;
		
			vector< vector<Data::P*> > node;
			
			inline vector<Data::P*> &operator()(int xi, int yi, int zi) {
				return node[xi + xn*(yi + yn*zi)];
			}
	
	} grid;
	
	
	Timer trf;
	double kfac;
	double r2_overlap;
	
	ICs(Data &data):
	data(data),
	kfac(0.977),
	r2_overlap(0.95*0.95)
	{ }
	
	void addHCPBox(int type1, int type2, Vector<double> R, Vector<int> D, Vector<double> S, Vector<double> V,  double t) {
			
		double dr = data.pw[type1][type2].s*kfac;
			
		for (int i=0; i<D.x; i++)
		for (int j=0; j<D.y; j++) 
			for (int k=0; k<D.z; k++){
					
				Data::P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				//pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
				pi.V *=0;
				
				//pi.R.x-= 0.5*S.x*dr;
				//pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
				//pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
				
				pi.R.x*= S.x;
				pi.R.y*= S.y;
				pi.R.z*= S.z;
				
				
				pi.R  += R;
				pi.V  += V;
				
				pi.type =  (i==0 || j==0 || k==0 || i==D.x-1 || j==D.y-1 || k==D.z-1) ? type2 : type1;

				data.p.push_back(pi);
					
					
			}
		data.p_size = data.p.size();
	}
	
	
	inline void find_overlaps(vector<Data::P*> &pi) {
	
	
		for (int i=0;i<pi.size()-1;i++) if (pi[i]->type!=-1)
		for (int j=i+1;j<pi.size();j++) {
						
			real dx = pi[i]->R.x - pi[j]->R.x; // 1
			real dy = pi[i]->R.y - pi[j]->R.y; // 1
			real dz = pi[i]->R.z - pi[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 < r2_overlap) {
				if (pi[i]->type==0) pi[j]->type=0;
				pi[i]->type =-1;
				continue;
			}
		}
	}


	inline void find_overlaps(vector<Data::P*> &pi,  vector<Data::P*> &pj) {
					
		for (int i=0;i<pi.size();i++) if (pi[i]->type!=-1)
		for (int j=0;j<pj.size();j++) {
						
			real dx = pi[i]->R.x - pj[j]->R.x; // 1
			real dy = pi[i]->R.y - pj[j]->R.y; // 1
			real dz = pi[i]->R.z - pj[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 < r2_overlap){
				if (pi[i]->type==0) pj[j]->type=0;
				pi[i]->type =-1;
				continue;
			}
				
			
		}
					
	}
	
	void remove_overlaps(){

		trf.set();
		
		grid.node.resize(grid.xn*grid.yn*grid.zn);
		
		for (int i=0; i<grid.node.size(); i++) grid.node[i].resize(0);
		
		
		double grid_05lenx = 0.5*grid.xn*grid.xs;
		double grid_05leny = 0.5*grid.yn*grid.ys;
		int k=0;
		data.p_size = data.p.size();
		for(int i=0; i<data.p_size; i++){
 			if (fabs(data.p[i].R.x) > grid_05lenx) { printf( "Warning grid_05lenx = %g < |%g = p[%i].R.x| we better stop!\n",i, grid_05lenx, data.p[i].R.x); exit(0); }
 			if (fabs(data.p[i].R.y) > grid_05leny) { printf( "Warning grid_05leny = %g < |%g = p[%i].R.y| we better stop!\n",i, grid_05leny, data.p[i].R.y); exit(0); }
 			
			grid(int((data.p[i].R.x+grid_05lenx)/grid.xs), int((data.p[i].R.y+grid_05leny)/grid.ys), int(data.p[i].R.z/grid.zs)).push_back(&data.p[i]);
			
		}

		
		 
		
		for (int zi=0; zi<grid.zn; zi++)
		for (int yi=0; yi<grid.yn; yi++)
		for (int xi=0; xi<grid.xn; xi++) {
			
			vector<Data::P*> &pi = grid(xi, yi, zi);
			
			if (pi.size()==0) continue;
							
			find_overlaps(pi);
			
			if (0 < xi && 0 < yi && zi+1<grid.zn)					find_overlaps(pi, grid(xi-1,yi-1,zi+1)); // 12
			if (0 < xi && zi+1<grid.zn)							   find_overlaps(pi, grid(xi-1,yi  ,zi+1)); // 1			
			if (0 < xi && yi+1<grid.yn)								find_overlaps(pi, grid(xi-1,yi+1,zi  )); // 13
			if (0 < xi && yi+1<grid.yn && zi+1<grid.zn)			find_overlaps(pi, grid(xi-1,yi+1,zi+1)); // 6
			
			if (0 < yi && zi+1<grid.zn)								find_overlaps(pi, grid(xi  ,yi-1,zi+1)); // 4
			if (zi+1<grid.zn)												find_overlaps(pi, grid(xi  ,yi  ,zi+1)); // 2
			if (yi+1<grid.yn && zi+1<grid.zn)						find_overlaps(pi, grid(xi  ,yi+1,zi+1)); // 5
			if (yi+1<grid.yn)												find_overlaps(pi, grid(xi  ,yi+1,zi  )); // 8
			
			if (xi+1<grid.xn && 0 < yi && zi+1<grid.zn)			find_overlaps(pi, grid(xi+1,yi-1,zi+1)); // 11
			if (xi+1<grid.xn && zi+1<grid.zn)						find_overlaps(pi, grid(xi+1,yi  ,zi+1)); // 3
			if (xi+1<grid.xn )											find_overlaps(pi, grid(xi+1,yi  ,zi  )); // 10
			if (xi+1<grid.xn && yi+1<grid.yn && zi+1<grid.zn)	find_overlaps(pi, grid(xi+1,yi+1,zi+1)); // 7
			if (xi+1<grid.xn && yi+1<grid.yn)						find_overlaps(pi, grid(xi+1,yi+1,zi  )); // 9
							
		}
	
		trf.get();
	
		for (int i=0; i<data.p_size; i++) 
			if (data.p[i].type == -1) data.p[i--] = data.p[--data.p_size];	
	
		grid.node.resize(0);

	}
	
	void save(string filename) {
			
		FILE *fp = fopen (filename.c_str() , "wb");
		
		struct MEM { Vector<float> R; int type;} *mem;
	
		mem = new MEM[data.p_size];
		for (int i=0; i<data.p_size; i++) {
			
			if (data.p[i].type!=0 && data.p[i].type!=1) printf("ICs::save() Warining!\n\tp[%i].type!=0 && p[%i].type!=1  = %i\n", i,i, data.p[i].type);
			
			mem[i].R = data.p[i].R;
			mem[i].type = data.p[i].type;
		}
		
		fwrite(mem, sizeof(MEM), data.p_size, fp);
		fclose(fp);
		delete [] mem;
	
	}
	
	void print(){
		
		printf("ICs::print()\n");
		
		printf("\tremoved particles: %i\n", data.p.size()-data.p_size);
		printf("\tfinding of overlaps: %g sec \n", trf.sin);
		printf("\tnumber of particle: %i\n", data.p_size);
		//printf("finding of overlaps: %g sec \n", trf.sin);
		//double g_memsize = ;
		printf("\tgrid.node[0]: %i byte\n", sizeof(grid.node[0]));
		printf("\tsizeof(grid.node[0]))*grid.node.size() = %g Mbyte\n", sizeof(grid.node[0])*(float)grid.node.size()/1024/1024);
		printf("\tgrid size: %g MByte\n", (sizeof(grid.node[0])*(float)grid.node.size() +  sizeof(Data::P*)*data.p.size())/1024/1024);
		printf("\n");
		//double p_memsize = sizeof(p[0])*double(p.size());
		
	}
	
	
};
